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The experimental control over the twist angle in twisted bilayer graphene has not been reported 
and its realistic structure is most likely incommensurate. In this paper, we develop a tight-binding 
virtual crystal approximation theory to study the electronic properties in incommensurate twisted 
bilayer graphene. The theory yields the electronic band structure and the local density of states 
for any incommensurate twist angle 9 between the graphene sheets. Angle dependent Van Hove 
singularities are observed in the numerically calculated local density of states. In accord with 
observations in scanning tunneling microscopy and spectroscopy, our theoretical calculation indicates 
that the rotation angle between graphene sheets does not result in a significant reduction in the Fermi 
velocity in comparison with monolayer graphene. The developed theory is quite general and can be 
applied to investigate the electronic properties in any incommensurate multilayer heterostructures. 
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The electronic structures of graphene bilayer systems 
depend considerably on their stacking structure W- 
The AB (Bernal) and AA stacking configurations, for 
example, exhibit very different electronic properties de¬ 
spite the subtle distinction between their stacking or¬ 
ders. More recently, the introduction of twisted bilayer 
graphene (tBLG), consisting of two sheets of graphene 
rotated relative to each other by an angle 9, has dramat¬ 
ically expanded the allotropes of the graphene bilayer 
and thin films. In particular, the tBLG demonstrated 
strong twist-dependent electronic spectra and properties 
[SHII], rendering the twist angle as a unique degree of 
freedom to modify the electronic properties of graphene 
multilayers. 

Twisted bilayer graphene is usually grown on the 
(0001) face of silicon-carbide (SiC) and its structure is 
theoretically classified as commensurate or incommen¬ 
surate, sensitively depending on the rotation angle [B]. 
Gommensurate tBLG is periodic and can be viewed as a 
honeycomb superlattice with more than four atoms per 
unit cell and larger lattice vectors. On the contrary, 
the incommensurate tBLG is non-periodic and charac¬ 
terized by structural disorder that changes with the in¬ 
commensurate twist angle. Commensuration in tBLG 
occurs theoretically at a discrete infinite set of twist an¬ 
gles m, with a vanishing probability that a randomly 
selected twist angle is commensurate. Experimentally 
fabricated tBLG are therefore most likely incommensu¬ 
rate, especially since the chemical synthesis of tBLG with 
a controlled rotation angle has not been reported unj. 

Despite the intensive theoretical effort [SHU, mm, 
a systematic theory to investigate the electronic struc- 
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ture in realistic incommensurate tBLG is still missing 
and highly indispensable. First-principles calculations in 
tBLG, for example, are restricted to commensurate tBLG 
|20) . [22]. Moreover, current tight-binding theories are 
restricted to the commensurate structures jSHTj, HMH], 
or, at their best, to incommensurate situations 
near commensuration [19] . 

In this Letter, we develop a non-trivial tight-binding 
virtual crystal approximation (TB-VCA) theory to study 
the electronic structure in incommensurate tBLG, repre¬ 
senting the first systematic study of the electronic prop¬ 
erties in this disordered system. The VGA theoretical 
approach is used in particular to determine the effective 
medium for a disordered tBLG system, derived in a direct 
manner from the quasi-infinite set of structural configura¬ 
tions for carbon atoms at the bilayer interface, for a given 
incommensurate twist angle. The theory hence yields 
the TB-VGA Hamiltonian matrix which turns out to be 
quasi-Hermitian as a consequence of the incommensurate 
twist angle 6. The hermiticity of the TB-VCA Hamilto¬ 
nian operator is successfully retrieved bytransposing the 
inner product, and eventually the Hilbert space associ¬ 
ated with the incommensurate tBLG system. The the¬ 
oretical TB-VCA model yields the electronic structure 
and the local electronic density of states (LDOS) for any 
incommensurate twist angle 9. 

Our TB-VCA theory underlines the pseudo-Hermitian 
representation of quantum mechanics p3H26] . Pseudo- 
Hermitian quantum mechanics is a recent and unconven¬ 
tional approach to quantum mechanics based on the use 
of non-Hermitian Hamiltonians, whose hermiticity can be 
retrieved by transposing the inner product defining the 
system Hilbert space. The pseudo-Hermitian representa¬ 
tion of quantum mechanics and the techniques developed 
in the course of its investigation have found applications 
in a variety of subjects [25], including condensed matter 
physics, scattering theory, relativistic quantum mechan¬ 
ics, and quantum cosmology. 

The TB-VCA theoretical method developed in this 
work for incommensurate graphene tBLG bilayer struc- 
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tures is quite general, and can be generalized in a 
straightforward manner to study the electronic structure 
and optoelectronic properties for any multiple incommen¬ 
surate twisted graphene heterostructure. 

The Letter is organized as follows. We present first the 
mathematical formalism developed for the TB-VCA the¬ 
ory. The numerical applications are then presented and 
discussed, and the last part is dedicated to conclusions 
and perspectives. 

The Tight-Binding (TB) representation is presented 
next for the particular tBLG problem. Consider the 
tBLG as composed of two graphene layers, denoted layer 
1 and 2, where layer 2 is rotated by an angle 9 with 
respect to layer 1. See Figjl]for the schematic represen¬ 
tation. Monolayer graphene is known to be composed 
of two sublattices named A and B. We hence refer to 
the sublattices in tBLG as Ai and Bi sublattices, with 
i = 1 and 2 for layers 1 and 2 respectively. Atoms in 
Ai and Bi sublattices are respectively labeled Ai{n) and 
Bi(n) with n = 0,l,2,.... In the absence of the periodic 
structure in incommensurate tBLG, our TB-VCA model 
considers referential coupled sites in each layer, denoted 
by {Ai(0), i3i(0)} and {^ 2 ( 0 ), 52(0)}, in layers 1 and 2 
respectively. 

In the VGA treatment of the disordered tBLG, the ef¬ 
fective medium is determined as a direct system average 
and the reciprocal-space techniques can hence be applied 
for the VGA periodic effective medium. Considering a 
single p 2 ;-orbital per site, the total electronic wave func¬ 
tion may hence be written as 

'f'k(i‘) = C'Ai'I'k;Ai (r) + C'Bj'I'k;Bi (r) 

+ + C'B2'I'k;B2(>^)> (1) 

where 4'k;a(i’) = [l/VV] - !■„), a = 

Ai,5i,^2,52, is the TB-VCA Bloch wave function for 
the a sublattice. Bold letters are used to denote vector 
quantities. JV is the number of unit cells in the crys¬ 
tal, Ca is the a-sublattice contribution to the total wave 
function, denotes the 2pz atomic orbital, Tq, is the po¬ 
sition of an a type atom with respect to a chosen origin, 
and k is the wave vector. 

We denote by V the vector space spanned by the basis 
{'I'k;Ai, 'I'k;Bi, ^'k;A2i 'f'k;B2 }. The standard inner prod¬ 
uct (.|.) defined on V is given by {x\4') = JfJ 
for any x and (j) in V. The vector space V endowed by the 
inner product (.|.) defines a Hilbert space T-L. This may be 
used to describe the incommensurate tBLG via the TB- 
VCA Hamiltonian operator H satisfying the Schrddinger 
equation 


iL4'k(r) = 5(k)4-k(r), (2) 

where 5(k) denotes the energy eigenvalue. It is impor¬ 
tant to note that the Hilbert space used to describe the 
incommensurate tBLG is not unique and may be changed 
by transposing the corresponding inner product. 

As in a standard tight-binding model, the energy eigen¬ 


value can be determined by solving the secular equation 

det [^"-^^(k)^] =0. (3) 

The label A denotes the energy bands. H and S 
denote respectively the Hamiltonian and the overlap 
matrices with matrix entities Sap = ('I'k;a|'I'k;/ 3 ) = 

///d3r4'^;^(r)Tk;/3(r) and Hap = (^^'k;a|-ff|^'k;/3) = 

///d3r4'*;^(r)ILTk;/3(r); a, l3 e {Ai, Bi, A 2 , B 2 }. We 
note that the definition of the matrices H and S follow 
directly the definition of the inner product in the Hilbert 
space. 

It is convenient to write each of the H and S in terms 
of 4 block matrices as follows 

[511(2x2) 7 Ji2(2x 2)1 
52i(2x2) 522(2x2) ’ 


where 


5^„(2 X 2) = 


Ha Ha mB„ 
Hb^a^ 5b„b„ 


( 5 ) 


and similarly for the overlap matrix S. 

The restriction to the single Pz - orbital per site yields 
ppn in-plane and ppa interlayer bonds in the incommen¬ 
surate tBLG system. In our calculations, we include up 
to the third nearest neighbors in-plane hopping and over¬ 
lap integrals, and the explicit form for the in-plane block 
matrices 5ii, H 22 , Su and S 22 can be found in refer¬ 
ence [27]. We further neglect the interlayer overlap in¬ 
tegrals, justified by the large perpendicular distance d± 
between the two graphene layers [dj_ « 3.35A). Gonse- 
quently S 12 = S 21 = 0. Note in this case that S*!! and 
S 22 are invariant under the twist angle. We are hence left 
with the determination of the interlayer block matrices 
5 i 2 and H 2 i- 

Through a numerical based statistical approach, we 
then analyzed the disordered interface morphology and 
the interlayer interactions in incommensurate tBLG. 
In particular, we have carefully investigated the inter¬ 
layer interactions between each of the referential sites 
{Ai(0), 5i(0)} from one layer interacting with sites 
{Aj{n),Bj{n)} in the other layer {j ^ i). 

Our study demonstrates the existence of a quasi¬ 
infinite set of configurations for carbon atoms at the 
tBLG bilayer interface, and hence a quasi-continuum set 
of interlayer bonds. This quasi-continuum set of inter¬ 
layer bonds requires a theoretical approach which yields 
the interlayer hopping integral as a function of the bond 
lengths. In our work, we choose the environment depen¬ 
dent tight-binding model, developed initially by Tang et 
al. [3S|, and optimized by Landgraf et al. [m, for the 
case of the tBLG. 

The key idea in our TB-VGA approach is hence to de¬ 
termine an ensemble average for the interlayer hopping 
integrals, and another average over the quasi-infinite en¬ 
semble of structural configurations. In this context, it 
is very useful to explore the possible statistical identi¬ 
ties over the structural disorder of the incommensurate 
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FIG. 1: The figure shows the numerically determined region 
Di in layer 1 (region D 2 in layer 2) which serve to simnlate 
the quasi-infinite ensembles of the strnctural configurations 
and interlayer bonds, at the incommensurate interface of the 
bilayer tBLG system, for different twist angles 6. This is done 
with respect to the referential sites ^i(O) (^12(0)) for the A 
sublattice. Similar figures may be schematically prepared for 
the other three configurations for the interactions of A and B 
sublattices in the two layers of the tBLG system. 


tBLG system. For example, our numerical analysis for 
the interlayer interactions demonstrates identical ensem¬ 
bles of interlayer bonds and configurations, whether we 
chose the referential sites ^i(O) or Bi(0) on layer 1 to in¬ 
teract with sites {^ 2 ( 71 ),i? 2 (tt)} in layer 2 or vice-versa. 
This analysis over the amorphous-like structural disorder 
of the incommensurate tBLG system, yields the following 
important identities 

HA\A2 ^A\B2 ^BiA2 ^BiB2: ( 6 ) 


and 


Note that the matrix HA 1 A 2 {HA 2 A 1 ) entities correspond 
respectively to the interlayer interactions between the ref¬ 
erential sites ^i(O) [^ 2 ( 0 )] with the A 2 [Ai] sublattices. 
This is also generically the case for the other identities 
of Eq|^ and Eqj^ 

The TB-VCA Hamiltonian matrix may hence be writ¬ 
ten from Eq. Eq. and Eq. in a reduced form 
as 


with 


_Hii HA 1 A 2 Q 
HA2A1Q H22 


( 8 ) 


Q 


I 1 
I 1 ’ 


(9) 


In what follows the TB virtual crystal approximation 
(TB-VGA) theory will be given. In order to proceed and 
determine the required TB-VGA ensemble averages, and 
consequently the matrix entities HaiA 2 and h is 

essential at this stage to simulate the ensembles of in¬ 
terlayer bonds and structural configurations. In Figj^ 
we present the numerically determined regions Di (in 
layer 1) and D 2 (in layers 2) which serve this purpose. 
For example, by varying the position of atom A 2 (0) such 
that its projection spans region Di, we formally simu¬ 
late the ensemble of interlayer bonds between the refer¬ 
ential site Ai(0) and sublattice A 2 for the ensemble of 
the structural configurations of the incommensurate in¬ 
terface. With this simulation, and using the environment 
dependent tight binding model approach, atom Ai(0) is 
found to interact with atoms {A 2 (n),n = 0, ...,6} of the 
A 2 sublattices. 

The matrix entity HaiA 2 may hence be determined in 
the TB-VCA approach as 

6 

HA 2 A 2 = '^i[M0)A2in)] 

n—0 

X exp [iirA2in) - r 2 i 2 (o)-k] , (10) 


HA2A1 — HA2B1 — Hb^Ai — Hb^Bi- 


(7) with 


t [A,(0)A,(0)] = 


flpjai^AAn) -rA.(0)||)r(j?niax- l|l~24,(ri) " fA. (0) 11) da:d^ 

max ||rA,(n) - I’yi.(o)ll) dxdy 


( 11 ) 


In Eqpr] II II designates the norm of a given vector, 
r„ is the position of atom a with respect to an ar¬ 
bitrary origin, F denotes the Heaviside step function. 


The i?max is the maximum interlayer separation be¬ 
yond which the interlayer interaction in neglected, t is 
the distance-dependent interlayer hopping integral de- 











4 



FIG. 2: The numerically calculated TB-VCA electronic structures for the tBLG system with incommensurate twist angles from 
left to right: 6 = 7° (A), 5° (B) and 2° (G), plotted along the high symmetry axes {FAi, KiMi, MiT} in the first BZ of layer 
1. Saddle points are observed in these electronic structures owing to the overlap between the Dirac cones. The circles represent 
the electronic band structure of a single graphene monolayer. The Fermi velocity near Ki for the tBLG system is observed to 
be identical to that of mass-less Dirac fermions in monolayer graphene. 


termined by the environment dependent tight binding 
model. t[Ai{ 0 )A 2 {n)] is the average interlayer hopping 
integral between atoms Ai(0) and A 2 {n). 

Due to symmetry (see FigQ, the averaged hopping 
integrals t[Ai{ 0 )A 2 {n)] for n = 1,...,6 are all identical. 
This yields hence 

Ha,a, = t[Aii0)A2m + 2t[Ai{0)A2il)] 

3 

X cos [(r^ 2 („) - rAi(o)) -k] . (12) 

n—1 

The same approach yields the matrix entity HA 2 A 1 as 
follows 


The quasi-Hermitian and Hermitian representations of 
the TB-VCA Hamiltonian operator H are formally equiv¬ 
alent [2S], [H], and yield the same results. The electronic 
structure for the incommensurate tBLG may hence be 
determined in a formal manner by solving the secular 
equation 

The TB-VCA numerical applications are presented 
next for different tBLG systems. The numerical values 
of the onsite energy, the in-plane hopping integrals, and 
the overlap integrals are adapted from reference [27] . The 
numerically calculated averaged interlayer hopping inte¬ 
grals are presented in Table 1. The numerical value of 
Rmax used in this calculation is i?max ~ 3.8A. 


HA 2 A, = t[A2(0)Ai(0)]+2t[A2(0)Ai(l)] 

3 

X Y cos [(rTi(n) - YA2{0)) -k] . (13) 

n—1 


Note that t [Ai(0)A2(0)] is invariant with twist angle, 
whereas t[Ai{ 0 )A 2 {n)] varies as an ensemble average 
with the twist angle. Also due to the bilayer inversion 
symmetry, t [Ai( 0 )A 2 (m)] = t[A 2 { 0 )Ai{m)] for m = 0, 
and m = 1, 2 ,.., 6 . 

By this, we complete the derivation of the Hamilto¬ 
nian matrix H, which turns out to be a quasi-hermitian 
matrix [25] . |26j . The hermiticity of the Hamiltonian 
operator H describing the incommensurate tBLG in the 
TB-VCA approach may be retrieved by transposing the 
inner product defined on the function vector space V. 

This is done by introducing a new inner product (.|.)^ 
on V, determined by a metric operator rj defined by 
fl'i’k-ai = _(1 + ^)5'k;ai and7)Tk;a2 = (l + ^k;a2 

with S = HA 2 A 1 /HA 1 A 2 Etnd a = A or B. This yields a 
Hermitian Hamiltonian matrix H' given as 


(1 -f 5)Hii (^HaiA2 + HA 2 A 1 ) Q 
{Ha,A2 + HA 2 A,) Q (1 + 1/^) H 22 

( 14 ) 


TABLE I: The ensemble averaged interlayer hopping integrals 
calculated numerically. 


f [Ai(0)A2(0)] (eV) f [Ai(0 )A2(1)] (eV) 


B 0 

1 ^ 

1 

0.2907 

0.1348 

6 = 5° 

0.2907 

0.1343 

to 

0 

0.2907 

0.1340 

e = r 

0.2907 

0.1339 


Our numerical calculations demonstrate the existence 
of two overlapped Dirac cones centered at the Dirac 
points Ki and K 2 in the Brillouin zones of layers 1 and 
2 respectively. In Figsj^ we plot the electronic struc¬ 
ture of incommensurate tBLG (solid red curves) along 
the high symmetry axes {FiFi, ATiMi, MiF} in the first 
BZ of layer 1. The figures from left to right correspond to 
the twist angles 0 = 7°, 5° and 2°. The open circles in the 
figures correspond to the electronic structure of graphene 
monolayer. The figures show that the Fermi velocity be¬ 
low the Van Hove singularities (VHS), is found in our 
theory to be identical to that of mass-less Dirac fermions 
in monolayer graphene. This result is in accordance with 
recent scanning tunneling microscopy and spectroscopy 
studies on incommensurate tBLG 0, m- 
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FIG. 3: The curves represent the numerically calculated 
LDOS in the TB-VCA model, in the neighborhood of the 
Dirac point Ki for the incommensurate twist angles 0 = 7°, 
5° and 2°, where the wider is the separation between the 
LDOS peaks the bigger is the twist angle. The tBLG LDOS 
present Van Hove singularities due to the overlap between the 
Dirac cones at Ki and K 2 - The low lying, almost flat, curve 
corresponds to the LDOS of single graphene monolayer. 

The electronic structures in Figsj^further demonstrate 
the existence of saddle points, induced from the band- 
overlap between the two Dirac cones. These saddle points 
generate the VHS peaks in the local density of states 
(LDOS) which have been observed experimentally [2], 
[Ms], [II]- The numerically calculated TB-VCA LDOS 


along the high symmetry axes {TKi, KiMi, MiT} are 
presented in Figs[^ for the twist angles 9 = 5°, 2° and 
1°. In accordance with experimental results, the saddle 
points and consequently the VHSs are observed to shift 
to lower energies when the twist angle is reduced. 

In conclusion we present the first systematic study of 
the electronic properties in incommensurate twisted bi¬ 
layer graphene. In particular, we have developed and 
applied the TB-VCA theory to determine the electronic 
band structure and local density of states in tBLG for 
any incommensurate twist angle 6. The theory estab¬ 
lishes for the disordered incommensurate tBLG system 
an averaged periodic VC A effective medium, deduced di¬ 
rectly from the quasi-continuum set of interlayer bounds 
and the quasi-infinite set of structural configurations 
characterizing the disordered system. The theory yields 
the quasi-Hermitian and Hermitian equivalent represen¬ 
tations for the TB-VCA Hamiltonian operator describ¬ 
ing the incommensurate tBLG. The numerical calcula¬ 
tion determines consequently the electronic band struc¬ 
ture and the angle dependent VHS in the LDOS for any 
incommensurate twist angle. The Fermi velocity is found 
to be identical to that of mass-less Dirac fermions in 
monolayer graphene. 

The current TB-VCA theory is quite general and can 
be applied to study the electronic structure in incom¬ 
mensurate twisted multilayer graphene and layered het¬ 
erostructures. Concerning perspectives, the authors are 
currently developing a dynamic variant of the non-local 
coherent potential approximation which fully describes 
the decay of electronic excitations in the tBLG system, 
adopting the current TB-VCA as a basic reference. 
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